Pattern formation in growing sandpiles 
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Adding grains at a single site on a flat substrate in the Abelian sandpile models produce beautiful 
complex patterns. We study in detail the pattern produced by adding grains on a two-dimensional 
square lattice with directed edges (each site has two arrows directed inward and two outward), 
starting with a periodic background with half the sites occupied. The size of the pattern formed 
scales with the number of grains added N as y/N. We give exact characterization of the asymptotic 
pattern, in terms of the position and shape of different features of the pattern. 



Many complicated and intricate patterns found in na- 
ture can be modelled by deterministic dynamics pQ. In 
Turing patterns [2] the final outcome is random due to 
the randomness in initial conditions. In the game of life 
[3] , one can get a very wide variety of patterns from sim- 
ple deterministic cellular automaton evolution rules, de- 
pending on the initial condition. 

While the real sand, poured at one point on a flat sub- 
strate produces a rather simple pyramidal shape, much 
more complex patterns are produced in the theoreti- 
cal models of sandpiles, like the Abelian sandpile model 
(ASM) [4]. Earlier studies have usually concentrated on 
determining the asymptotic shape of the growing clus- 
ter [5j[6]. Other special configurations in the model, like 
the identity [7] , or the stable state produced from special 
unstable states also show complex internal self-similar 
structures [8]. The limiting shape has been determined 
in the related rotor-router model, and the model of di- 
visible sandpiles with multiple sites of addition [9]. 

In this paper, we study the asymptotic pattern pro- 
duced by adding N grains of sand at a single site on a 
two dimensional Abelian sandpile model starting from a 
periodic background, and allowing the system to relax. 
It is easy to see that the diameter of the pattern grows 
as y/N. Interestingly, for large TV, the pattern shows a 
proportionate growth, with different parts of the pattern 
all growing as y/N. This is thus different from earlier- 
studied models of growth such as diffusion limited aggre- 
gation, Eden model etc. [10 , where the growth occurs 
mainly at the surface. 

The standard square lattice produces a rather compli- 
cated pattern (FigfTJi), and it has not been possible to 
characterize it so far. We consider two variations, assign- 
ing orientations to the edges of the lattice, as shown in 
FigJSJi and ^). The initial state was chosen to be a pe- 
riodic checkerboard arrangement of sites with heights 
and 1. The asymptotic pattern produced in the two cases 
turns out to be the same, and is shown in Fig [if). Taking 
some qualitative features of the observed pattern ( e.g. 
only two types of patches are present, and they are all 
3- or 4- sided polygons) as input, we show how one can 
get a complete and quantitative characterization of the 
pattern. We show that the pattern has exact 8-fold ro- 
tational symmetry, and determine the exact coordinates 





FIG. 1: Stable configurations for the Abelian sandpile model 
obtained by adding particles at one site, (a) Undirected 
square lattice, initial configurations with all heights 2, and 
2 x 10 5 particles added, color code: red=0, blue=l, green=2, 
yellow=3. (b) The F-lattice of Fig]^ with initial checker- 
board configuration, with 2 x 10 5 particles added, color code: 
green=0, yellow=l. The apparent green regions in the picture 
represent the patches with checkerboard configuration. 



FIG. 2: The directed square lattices studied in this paper (a) 
the F-lattice (b) the Manhattan lattice. 



of all the boundaries in the asymptotic pattern. We dis- 
cuss some other cases, where exactly the same pattern is 
obtained. 

In the two lattices we studied (Figj2|, each bond of 
the lattice is directed with two in- arrows, and two out- 
arrows at each vertex. The ASM on these is defined by 
the toppling rule: A site (x, y) is unstable if the number 
of grains at the site z x ^ y > 2, and then transfers one grain 
each in the direction of its outward arrows. We start with 
an initial configuration in which z XjV = 1, for sites with 
(x + y) = even, and otherwise. 
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We used a lattice large enough so that no avalanches 
started from the origin reach the boundary. Using the 
Abelian property, we add all N particles in the beginning, 
and relax the configuration to get the final pattern. The 
result of adding N = 2 x 10 5 particles on the F-lattice is 
shown in FigjlJ). The pattern formed on the Manhattan 
lattice is indistinguishable at large scales. The pattern 
is identical to Fig[TJ), except that the thin lines of l's 
forming two triangles outside the octagon are rotated by 
45° in the Manhattan case. Since the lattices are quite 
different, this is quite intriguing. 

We start by setting up some general theoretical frame- 
work, which is independent of the details of the partic- 
ular lattices studied. Formally, we can characterize the 
asymptotic pattern in terms of the rescaled coordinates, 
£ = x/\/N, tj = y/VN and the density function 77) 
which gives the local density of grains in the pattern in 
a small rectangle of size A£, At? about the point (£,ry), 
with TV" 1 / 2 <C Af , Ar; <C 1. 

Equivalently, we can describe the asymptotic pattern 
in term of the rescaled toppling function (/>(£, 77). Let 
Tjy(x,y) be the numbers of toppling at site (x,y) when 
N particles are added at the origin, and the configuration 
is relaxed. We define 

M, V ) = lim ±T N (lVNZ\, [Vn V \) (1) 

N^oo ZiV 

where floor function [x\ is the largest integer less than 
x. From the conservation of sand grains, it is easily seen 
that 77) is related to the density function p(£, 77) by 

(|* + l^' 7 ?) = - 6®5( V ) (2) 

where excess density Ap(£, 77) is the difference between 
p(£, 77) and the initial density po(^T]). 

It was already noted [IT] that for N large p(£, 77) tends 
to a nontrivial limit, and the asymptotic pattern is made 
of distinct regions, called 'patches'. Typically inside a 
patch the heights are periodic in space, and there are 
few defect-lines, which move with iV, but do not change 
the macroscopic density p(£, 77). Then, the coarse grained 
function p(£,ry) takes constant rational value in each 
patch. Also in each patch of constant Ap(£, 77), 77) is 
a quadratic function, and was first noted in [TT] , We in- 
dicate the proof here. For all patches the function 77) 
is Taylor expandable around any point inside the patch. 
Consider any term of order > 3 in the expansion, for ex- 
ample the term ~ (A£) 3 . This can only arise due to a 
term ~ (Ax) 3 /VN in T(x,y). Then the integer func- 
tion T(x,y) will change discontinuously at intervals of 
Ax ~ 0(N 1 ^ 6 ) leading to infinitely many defect-lines 
in the asymptotic pattern. However there are no such 
feature in Fig(T]2 or Fig[TJ). Therefore inside a patch of 
constant Ap(£, 77), 77) can at most be quadratic in 
£ and 77, and in each periodic patch, the toppling func- 
tion T(x,y) is sum of two terms: a part that is a sim- 
ple quadratic function of x and and a periodic part. 



The periodic part averages to zero, and does not con- 
tribute to the coarse-grained function (/>(£, 77). In some 
patterns, there are regions of finite fractional area which 
show aperiodic height patterns. In these regions </>(£, 77) 
is not quadratic and are harder to characterize. 

Now consider two neighboring periodic patches P and 
P' with mean densities p and p' respectively. Let the 
quadratic toppling function be Q(£, 77) and Q f (^r]) in 
these patches. Then the boundary between the patches 
is given by the equation 77) = Q f (^r]). As the 
derivatives of (j) are also continuous across the bound- 
ary, the boundary between two periodic patches must be 
a straight line, and 

Q'(z,v) = Q(t;,v) + \(p , -p)ii (3) 

where l± is the perpendicular distance of (£, 77) from the 
boundary. We can start with a periodic patch P, and 
go to another patch P' by more than one path. Since 
the final quadratic function at P' should be the same 
whichever path we take, this imposes consistency con- 
ditions which restricts the allowed values of slopes of 
boundaries. Consider a point zo where n periodic patches 
meet, with n > 2 (Fig {3)2). If the jth boundary at this 
point makes an angle Oj with the x-axis, and the density 
of the patch in the wedge Oj < < Oj+i is (Figj3ji) 
then using Eqj3] repeatedly for all n patches around zq 
we get that the following equation must be true for all 0: 

n 

^2(p^ 1 -p J )sin 2 (0-0 J ) = 1 (4) 
with p n +i = p n . This is equivalent to the condition: 

n 

^2(p J+1 -p 3 )e^=0 (5) 
3=1 

For n = 3, with pi 7^ p2 7^ P3, this equation has only 
trivial solutions with Oj equal to or tt for all j. Hence, 
only n > 4 are allowed. 

We now discuss how the exact function 77) can be 
determined for our problem. We note that in PigJTJ>, 
there are no aperiodic patches, only two types of periodic 
patches, where 77) only take values 1 or 1/2. Also, 
the slopes of the boundaries between patches only take 
values 0, ±1, 00. The patches are typically dart shaped 
quadrilaterals, and some triangles (which may be consid- 
ered as degenerate quadrilaterals with one side of length 
zero). These simplifications, not present in FigjTJi, make 
possible a full characterization of the pattern in Fig|TJ). 

Given that there are only these two types of patches, 
we only need to look for possible patterns where Ap takes 
piece wise constant values 1/2 or 0. From Eq.([2|, we see 
that we can think of 77) as the potential produced 
by a point charge at the origin, and a charge cloud with 



FIG. 3: (a) n different periodic patches of density pi,...,p n 
meeting at point zq. (b) The pattern in FigJTJ? is obtainable 
by putting together square tiles of different sizes. Each of the 
tiles is divided into two halves of different density. 

areal density — Ap(£,7/), with total charge zero. The ba- 
sic principle which selects the actual stable pattern out 
of many is a version of the principle of minimum dissi- 
pation: It is a stable state reached by minimum number 
of toppling. (This follows immediately from the toppling 
rules, where no toppling occurs unless forced). 

The requirement that 77) be exactly zero, in the 
region outside the pattern, implies that all the multi- 
pole moments of the charge distribution Ap(£, 77) are 
exactly zero. We show below that the conditions that 
Ap takes only two values, the potential function is ex- 
actly quadratic within a periodic patch, and the slopes 
of the boundaries are only 0, ±1, 00, fix the allowed pat- 
tern uniquely. 

We start by determining the exact asymptotic size of 
the pattern. We note from Fig{T}> that the boundary of 
the pattern is an octagon ( we shall prove later that this 
is a regular octagon ). In fact there are four lines of l's 
outside the octagon. But these has zero areal density in 
the limit N —> 00, and do not contribute to 77). We 
will ignore these in the following discussion. 

Let B be the minimum boundary square containing 
all (£, 77) that have a non-zero charge density 77). We 
observe that B can be considered as a union of disjoint 
smaller squares, each of which is divided by diagonal into 
two parts where Ap(£, 77) takes values 1/2 and [Figj3f)]. 
This is seen to be true for the outer layer patches. To- 
wards the center, the squares are not so well resolved. 
Assuming that this construction remains true all the way 
to the center, in the limit of large TV, the mean density 
of the negative charge in the bounding square = 1/4. 
Given that the total amount of negative charge is — 1, 
the area of the bounding square should be 4. Hence we 
conclude that the equation of the boundary of the mini- 
mum bounding square are 

ICI = 1, H = 1 (6) 

Let Nb be the minimum number of particles that have to 
be added so that at least one site at y = b topples. We 
find that for b = 10, 50, 100, and 300, ^/W h = 10.770, 
49.436, 98.894 and 297.798. Clearly the boundary dis- 
tance b tends to y/N for large N. 



FIG. 4: Two representations of the adjacency graph of the 
pattern. Here the vertices are the patches, and the edges 
connect the adjacent patches, (a) Representation as a planar 
graph (b) as a graph of wedge of angle 4-zr formed by glueing 
together the eight quadrant graphs at the origin. 



We now describe the topological structure of the pat- 
tern. This is characterized by its adjacency graph 
[Fig j4jx] , where each vertex denotes a patch, and a bond 
between the vertices is drawn if the vertices share a com- 
mon boundary. It is convenient to think of the triangular 
patches in the pattern as degenerate quadrilaterals, with 
one side of length zero. Then we see that the adjacency 
graph is planar with each vertex of degree four, except a 
single vertex of coordination number eight corresponding 
to the exterior of the pattern. The graph has the struc- 
ture of a square lattice wedge, with wedge angle 47r. The 
square lattice structure of the adjacency graph is seen 
most directly by applying a^' = l/z 2 transformation to 
the picture (used earlier in [H]), where z = £ + irj, and 
view it in the complex z'-plane. Thus, one can equiva- 
lently represent the graph as a square grid on a Riemann 
surface of two sheets (figj4f>). 

We now use the qualitative information obtained from 
the adjacency matrix of the observed pattern, to obtain 
quantitative prediction of the exact coordinates of all the 
patches. Consider an arbitrary patch P, having an excess 
density 1/2. The potential function in the patch is a 
quadratic function of (£, 77) and we parametrize it as 

+<U + e P 77 + / P (7) 

The potential function in a patch P having zero excess 
density will be parametrized as 

0p(&*7) = ^p(^ 2 -^ 2 ) + ^P^ + ^P^ + e p 77 + / p (8) 

Now consider two neighboring patches P and P' with 
excess densities 1/2 and respectively. Then using the 
matching condition Eq.(3), it is easy to show that if the 
boundary between them is a horizontal line 77 = 77 p , we 
must have 

77V = m r + 1; n p /=n p , dp/ = d p 
e p / = e p +77 P /2, f p ,=f p -r£/4 (9) 
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Similar calculation for other boundaries show that across 
a vertical boundary, going from a patch of higher density 
to the one of lower density, we have Am p = — 1 and 
An p = 0. Across a boundary with slope ±1, Am p = 0, 
and An p = ±1. 

In the outermost patch, clearly rj) = 0, and for 
this patch both m and n are zero. It follows that all ra p 
and n p take integer values. In the following, we denote 
a patch by integers (m,n), and write the corresponding 
coefficients <i p , e p , and f p as d m>n , e m?n and f m ,n- With 
this convention, the matching conditions in Eq.Q can 
be rewritten as 

4+i,n <^m,n, e m+ i )n e m)n 7?m,n/2, (m + n) odd 

(10) 

Using similar matching conditions for the boundary of 
patch (m, n) with slope ±1, we get the conditions 

dm,n+i ~ dm, n = e m , n - e m , n+ i , (m + n) odd 
d m ,n-i - ^m,n = e m , n _i - e m , n , (m + n) odd (11) 

We can eliminate the variables <i m , n and e m5n with (m + 
n) even using Eq.(10) and Eq.(ll). Then the equations 
become 

^m+2,n ^m,n ?7ra,n/2 (12) 

dm—2,n ^m,n — £ra,n/2 (13) 
dm—l,n—l d mn 6 m _|_i ?n _i 6rn,n (14) 

dm— l,n+l d m ^ n — [^m+l ^m,n] (15) 

It is convenient to introduce the complex variables z = 
£ + Z77, M = m + in and D = d + ie. In these variables 
we can write Eq.(7) as 



|zz+|ite[z 2 M + £>*]+/ (16) 

O O 



Under a rotation of axes by an angle 0, z — > z' = ze z6> , 
the requirement that is invariant is satisfied if we have 



M' = Me 2ie : D' = De l 



(17) 



On the (m, n) lattice, with (m + n) odd, the natural 
basis vectors are (1,1) and (1,-1). Let us call these a 
and (3. We define the finite difference operators A± a and 
A ±/ 3 by 

A ±a f(z) = f(z±a)-f(z) 

A ±(3 f(z) = f(z±(3)-f(z) (18) 

Then the equations (14-15) can be written as 

A- a d = Ape 
A_f3d=-A a e (19) 

These equations are the discrete analog of the famil- 
iar Cauchy-Riemann conditions connecting the partial 
derivatives of real and imaginary parts of an analytic 



function where the role of the analytic function is played 
by D — d + ie. 

From Eq.(14) and Eq.(15), it is easy to deduce that D 
satisfies the discrete Laplace's equation 



[A a A_ a + A p A_p]D = 



(20) 



If m and n are large, the corresponding patch is near 
the origin (|£| + |^| is small), and where the leading be- 
havior of 0(£,r/) is given by </>(£, 7?) ~ -^log(£ 2 + rj 2 ). 
Consider a point Zq, such that at zq 

d 2 4>/d£ 2 * m/4; d 2 4>/d^dn w n/4, (21) 

Then, zo would be expected to lie in the patch labeled 
by (m,n). This gives z « ±(ttM/2)" 1/2 . Then, setting 
d(j)/dz equal to 5/2 gives us 



1 



2tt 



(22) 



The equation (20), subjected to the behavior at large 
\n\ given by Eq.(22) on the 4tt- wedge graph (for 



m 



each value of (m, n), D m ^ n has two values) has an unique 
solution. Clearly the solution has eight fold rotational 
symmetry about the origin in the (m, n) space. This 
implies that 



,•1/2 



L> m)U ; for all (m,n). 



(23) 



Given D m ^ n , its real and imaginary parts determine d m ^ n 
and e m?n , and using Eq.(12, 13) we determine the exact 



positions of all the patch corners. The exact eight-fold ro- 
tational symmetry of the adjacency graph of the pattern, 
and the fact that D satisfies Eq.(20) on the adjacency 
graph together imply the eight-fold rotational symmetry 
of all the distances in the pattern. 

We have not been able to find a closed-form formula 
for D m ^ n . But the system of coupled linear equations 



(20) can be determined numerically to very good pre- 
cision by solving it on a finite grid — L < m, n < L 



with the condition in Eq.(22) imposed exactly at the 
boundary. We determined d m?n and e m:n numerically 
for L = 100, 200, 400, and extrapolated our results for 
L -> oo. We find d h0 = 0.5000 and d 2 ,i = 0.6464, in 
perfect agreement with the exact theoretical values 1/2 
and 1 — l/2>/2 respectively. 

Our arguments above can be extended to other two 
dimensional lattices, so long as there are only two allowed 
values of Ap. While this is not clear why, this seems 
to happen for the Manhattan lattice (Fig|2^), for initial 
density 1/2. Also, this happens on the F-lattice, with 
a periodic background pattern with initial density 5/8 
[zij = 1 if i + j even, or congruent to (0, 1) or 

(2, 3) mod 4]. In some other cases, like the F-lattice, 
with initially all sites empty, the pattern is very similar, 
but there are some non periodic patches in the outermost 
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ring. Since the behavior of rf) in such patches is not 
known, the equations for D m ^ n do not close in this case. 
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